##### ####################################################
#####                                               ######
#####           Responsiveness Analysis            
#####                                               ######
##### ####################################################

# init ------------------------------------------------------------

rm(list=ls())
set.seed(221186)

# Load libraries

library(data.table) # 1.11.4
library(mgcv) # 1.8.24
library(stargazer) # 5.2.2
library(xtable) # 1.8.2

source("utils.R")

# Load data

load("working/speeches_influence.Rdata")

### ################################################
### Prepare data
### ################################################

speeches <- speeches[minister_present == T]

### ################################################
### Models
### ################################################

## OLS -- interaction specification

speeches$yearmon.numeric <- as.numeric(speeches$yearmon)-1997

mod1 <- as.formula("sim.words ~ Gender*minister_gender + next.speech.same.party")
mod2 <- as.formula("sim.words ~ Gender*minister_gender + next.speech.same.party + debate_department")
mod3 <- as.formula("sim.words ~ Gender*minister_gender + next.speech.same.party + as.factor(yearmon)")
mod4 <- as.formula("sim.words ~ Gender*minister_gender + next.speech.same.party + debate_department + as.factor(yearmon)")
mod5 <- as.formula("sim.words ~ Gender*minister_gender + next.speech.same.party + debate_department*yearmon.numeric + as.factor(yearmon)")
mod6 <- as.formula("sim.words ~ Gender*minister_gender + next.speech.same.party + debate_department*poly(yearmon.numeric,2) + as.factor(yearmon)")

model.data <- speeches[minister_in_debate==F & next.speech.minister==T ]

lm.mod1 <- cluster.function(mod1, cluster = "debate_department", data = model.data)
lm.mod2 <- cluster.function(mod2, cluster = "debate_department", data = model.data)
lm.mod3 <- cluster.function(mod3, cluster = "debate_department", data = model.data)
lm.mod4 <- cluster.function(mod4, cluster = "debate_department", data = model.data)
lm.mod5 <- cluster.function(mod5, cluster = "debate_department", data = model.data)
lm.mod6 <- cluster.function(mod6, cluster = "debate_department", data = model.data)

model.data$debate_department_factor <- as.factor(model.data$debate_department)
gam.mod.out <- gam(sim.words ~ Gender*minister_gender + next.speech.same.party + s(yearmon.numeric, by=debate_department_factor) + debate_department_factor + as.factor(yearmon), data=model.data)
names(gam.mod.out$coefficients)[1] <- "Constant"
names(gam.mod.out$coefficients)[grep("GenderF:minister_genderF",names(gam.mod.out$coefficients))] <- "Interaction"
gam.se <- summary(gam.mod.out)$se
names(gam.se)[1] <- "Constant"
names(gam.se)[grep("GenderF:minister_genderF",names(gam.se))] <- "Interaction"

party.fe <- c("Same Party Control","$\\checkmark$","\\checkmark","$\\checkmark$","\\checkmark","\\checkmark","\\checkmark","\\checkmark")
ministry.fe <- c("Ministry FEs","$\\times$","\\checkmark","$\\times$","\\checkmark","\\checkmark","\\checkmark","\\checkmark")
month.fe <- c("Month FEs","$\\times$","$\\times$","\\checkmark","\\checkmark","\\checkmark","\\checkmark","\\checkmark")
linear.tt <- c("Linear time trends","$\\times$","$\\times$","$\\times$","$\\times$","\\checkmark","\\checkmark","$\\times$")
quadratic.tt <- c("Quadratic time trends","$\\times$","$\\times$","$\\times$","$\\times$","$\\times$","\\checkmark","$\\times$")
gam.tt <- c("Flexible time trends","$\\times$","$\\times$","$\\times$","$\\times$","$\\times$","$\\times$","\\checkmark")


sink("latex/tables/responsiveness_tables.tex")
mod_stargazer(lm.mod1, lm.mod2, lm.mod3, lm.mod4, lm.mod5, lm.mod6, gam.mod.out,
              keep=c("Gender","minister_gender","Interaction","Constant"), 
              order=c("Gender","minister_gender","Interaction","Constant"), 
              se=list(lm.mod1$clustered.se, lm.mod2$clustered.se, lm.mod3$clustered.se, lm.mod4$clustered.se, lm.mod5$clustered.se, lm.mod6$clustered.se,gam.se), 
              covariate.labels=c("Female","Female Minister","Interaction","Constant"), 
              keep.stat=c("n","rsq","adj.rsq"), 
              dep.var.caption="", 
              dep.var.labels="\\emph{Responsiveness}", 
              add.lines = list(party.fe, ministry.fe, month.fe, linear.tt, quadratic.tt,gam.tt),
              column.sep.width="0.25pt",
              label="tab:ols_cosine", 
              title="Responsiveness (OLS)", 
              no.space=T, model.names=F)
sink()

## F-tests

mod1f <- as.formula("sim.words ~ next.speech.same.party")
mod2f <- as.formula("sim.words ~ next.speech.same.party + debate_department")
mod3f <- as.formula("sim.words ~ next.speech.same.party + as.factor(yearmon)")
mod4f <- as.formula("sim.words ~ next.speech.same.party + debate_department + as.factor(yearmon)")
mod5f <- as.formula("sim.words ~ next.speech.same.party + debate_department*yearmon.numeric + as.factor(yearmon)")
mod6f <- as.formula("sim.words ~ next.speech.same.party + debate_department*poly(yearmon.numeric,2) + as.factor(yearmon)")

lm.mod1f <- cluster.function(mod1f, cluster = "debate_department", data = model.data)
lm.mod2f <- cluster.function(mod2f, cluster = "debate_department", data = model.data)
lm.mod3f <- cluster.function(mod3f, cluster = "debate_department", data = model.data)
lm.mod4f <- cluster.function(mod4f, cluster = "debate_department", data = model.data)
lm.mod5f <- cluster.function(mod5f, cluster = "debate_department", data = model.data)
lm.mod6f <- cluster.function(mod6f, cluster = "debate_department", data = model.data)

f_stats <- c(anova(lm.mod1f, lm.mod1)$F[2],
  anova(lm.mod2f, lm.mod2)$F[2],
  anova(lm.mod3f, lm.mod3)$F[2],
  anova(lm.mod4f, lm.mod4)$F[2],
  anova(lm.mod5f, lm.mod5)$F[2],
  anova(lm.mod6f, lm.mod6)$F[2])


## OLS -- split sample specification

mod1_split <- as.formula("sim.words ~ minister_gender + next.speech.same.party")
mod2_split <- as.formula("sim.words ~ minister_gender + next.speech.same.party + debate_department")
mod3_split <- as.formula("sim.words ~ minister_gender + next.speech.same.party + as.factor(yearmon)")
mod4_split <- as.formula("sim.words ~ minister_gender + next.speech.same.party + debate_department + as.factor(yearmon)")
mod5_split <- as.formula("sim.words ~ minister_gender + next.speech.same.party + debate_department*yearmon.numeric + as.factor(yearmon)")
mod6_split <- as.formula("sim.words ~ minister_gender + next.speech.same.party + debate_department*poly(yearmon.numeric,2) + as.factor(yearmon)")
mod7_split <- as.formula("sim.words ~ minister_gender + next.speech.same.party + s(yearmon.numeric, by=debate_department_factor) + debate_department_factor + as.factor(yearmon)")

model.data.male <- speeches[minister_in_debate==F & next.speech.minister==T & Gender == "M"]
model.data.male$debate_department_factor <- as.factor(model.data.male$debate_department)
model.data.female <- speeches[minister_in_debate==F & next.speech.minister==T & Gender == "F"]
model.data.female$debate_department_factor <- as.factor(model.data.female$debate_department)

lm.mod1.male <- cluster.function(mod1_split, cluster = "debate_department", data = model.data.male)
lm.mod2.male <- cluster.function(mod2_split, cluster = "debate_department", data = model.data.male)
lm.mod3.male <- cluster.function(mod3_split, cluster = "debate_department", data = model.data.male)
lm.mod4.male <- cluster.function(mod4_split, cluster = "debate_department", data = model.data.male)
lm.mod5.male <- cluster.function(mod5_split, cluster = "debate_department", data = model.data.male)
lm.mod6.male <- cluster.function(mod6_split, cluster = "debate_department", data = model.data.male)
gam.mod.out.male <- gam(mod7_split, data = model.data.male)
names(gam.mod.out.male$coefficients)[1] <- "Constant"
gam.male.se <- summary(gam.mod.out.male)$se
names(gam.male.se)[1] <- "Constant"

sink("latex/tables/responsiveness_tables_male.tex")
mod_stargazer(lm.mod1.male, lm.mod2.male, lm.mod3.male, lm.mod4.male, lm.mod5.male, lm.mod6.male, gam.mod.out.male,
              keep=c("minister_gender","Constant"), 
              order=c(2,1), 
              se = list(lm.mod1.male$clustered.se, lm.mod2.male$clustered.se, lm.mod3.male$clustered.se, lm.mod4.male$clustered.se, lm.mod5.male$clustered.se, lm.mod6.male$clustered.se, gam.male.se), 
              covariate.labels=c("Female Minister","Constant"), 
              keep.stat=c("n","rsq","adj.rsq"), 
              dep.var.caption="", 
              dep.var.labels="Responsiveness (male MPs)", 
              add.lines = list(party.fe, ministry.fe, month.fe, linear.tt, quadratic.tt,gam.tt),
              column.sep.width="0.25pt",
              label="tab:ols_cosine", 
              title="Responsiveness (OLS)", 
              no.space=T, model.names=F)
sink()


lm.mod1.female <- cluster.function(mod1_split, cluster = "debate_department", data = model.data.female)
lm.mod2.female <- cluster.function(mod2_split, cluster = "debate_department", data = model.data.female)
lm.mod3.female <- cluster.function(mod3_split, cluster = "debate_department", data = model.data.female)
lm.mod4.female <- cluster.function(mod4_split, cluster = "debate_department", data = model.data.female)
lm.mod5.female <- cluster.function(mod5_split, cluster = "debate_department", data = model.data.female)
lm.mod6.female <- cluster.function(mod6_split, cluster = "debate_department", data = model.data.female)
gam.mod.out.female <- gam(mod7_split, data = model.data.female)
names(gam.mod.out.female$coefficients)[1] <- "Constant"
gam.female.se <- summary(gam.mod.out.female)$se
names(gam.female.se)[1] <- "Constant"

sink("latex/tables/responsiveness_tables_female.tex")
mod_stargazer(lm.mod1.female, lm.mod2.female, lm.mod3.female, lm.mod4.female, lm.mod5.female, lm.mod6.female, gam.mod.out.female,
              keep=c("minister_gender","Constant"), 
              order=c(2,1), 
              se=list(lm.mod1.female$clustered.se, lm.mod2.female$clustered.se, lm.mod3.female$clustered.se, lm.mod4.female$clustered.se, lm.mod5.female$clustered.se, lm.mod6.female$clustered.se, gam.female.se), 
              covariate.labels=c("Female Minister","Constant"), 
              keep.stat=c("n","rsq","adj.rsq"), 
              dep.var.caption="", 
              dep.var.labels="Responsiveness (female MPs)", 
              add.lines = list(party.fe, ministry.fe, month.fe, linear.tt, quadratic.tt,gam.tt),
              column.sep.width="0.25pt",
              label="tab:ols_cosine", 
              title="Responsiveness (OLS)", 
              no.space=T, model.names=F)
sink()


## ####################################
## Plot effect sizes
## ####################################

sum.coefficients.ci <- function(model, mod.vcov=NULL, var2 = "^minister_genderF$", var1 = "^GenderF", int = "^Interaction$", baseline.var2, baseline.int){
  
  if(is.null(mod.vcov)) mod.vcov <- vcov(model)
  
  mod.coef <- na.omit(coef(model))
  colnames(mod.vcov) <- names(mod.coef)
  rownames(mod.vcov) <- names(mod.coef)
  
  ## Find the variables in the covariance matrix
  var1.position <- grep(var1, colnames(mod.vcov))
  var2.position <- grep(var2, colnames(mod.vcov))
  interaction.position <- grep(int, colnames(mod.vcov))
  
  # Find the variances
  var1.var <- diag(mod.vcov)[var1.position]
  var2.var <- diag(mod.vcov)[var2.position]
  interaction.var <- diag(mod.vcov)[interaction.position]
  
  # Find the covariance
  var2var3.covariance <- mod.vcov[var2.position,interaction.position]
  
  # Point esitmates for b2 and b3
  var1.pe <- mod.coef[var1.position]
  var2.pe <- mod.coef[var2.position]
  interaction.pe <- mod.coef[interaction.position]	
  
  # Sum the point estimates for b2 and b3
  pe.for.sum <- mod.coef[var2.position] + mod.coef[interaction.position]
  
  # Find the standart error
  se.for.sum <- sqrt(var2.var + interaction.var + 2*(var2var3.covariance))
  
  # Find the confidence intervals
  var1.upper <- var1.pe + 1.96 * sqrt(var1.var)
  var1.lower <- var1.pe - 1.96 * sqrt(var1.var)
  var2.upper <- var2.pe + 1.96 * sqrt(var2.var)
  var2.lower <- var2.pe - 1.96 * sqrt(var2.var)
  interaction.upper <- interaction.pe + 1.96 * sqrt(interaction.var)
  interaction.lower <- interaction.pe - 1.96 * sqrt(interaction.var)
  sum.upper <- pe.for.sum + 1.96 * se.for.sum
  sum.lower <- pe.for.sum - 1.96 * se.for.sum
  
  # Is the interaction significant?
  sig.test  <- coeftest(model, vcov = mod.vcov)
  sig.test <- sig.test[interaction.position,4]<0.05
  
  effect.size.var2 <- var2.pe/baseline.var2
  effect.size.var2.lower <- var2.lower/baseline.var2
  effect.size.var2.upper <- var2.upper/baseline.var2
  
  effect.size.int <- pe.for.sum/baseline.int
  effect.size.int.lower <- sum.lower/baseline.int
  effect.size.int.upper <- sum.upper/baseline.int
  
  # Bind it all together
  out.coef <- data.frame(var1.pe, var1.se= sqrt(var1.var), var1.lower, var1.upper,var2.pe, var2.se = sqrt(var2.var), var2.lower, var2.upper, pe.for.sum, se.for.sum, sum.lower, sum.upper, interaction.pe, interaction.se = sqrt(interaction.var), interaction.upper, interaction.lower, effect.size.var2 , effect.size.var2.lower, effect.size.var2.upper, effect.size.int , effect.size.int.lower, effect.size.int.upper, sig.test)
  return(out.coef)
}

baseline.male <- model.data[minister_in_debate==F & minister_gender=="M" & Gender=="M",mean(sim.words,na.rm=T)]
baseline.female <- model.data[minister_in_debate==F & minister_gender=="M" & Gender=="F",mean(sim.words,na.rm=T)]

model1est <- sum.coefficients.ci(model = lm.mod1, mod.vcov = lm.mod1$clus.vcov, baseline.var2 = baseline.male, baseline.int = baseline.female)
model2est <- sum.coefficients.ci(model = lm.mod2, mod.vcov = lm.mod2$clus.vcov, baseline.var2 = baseline.male, baseline.int = baseline.female)
model3est <- sum.coefficients.ci(model = lm.mod3, mod.vcov = lm.mod3$clus.vcov, baseline.var2 = baseline.male, baseline.int = baseline.female)
model4est <- sum.coefficients.ci(model = lm.mod4, mod.vcov = lm.mod4$clus.vcov, baseline.var2 = baseline.male, baseline.int = baseline.female)
model5est <- sum.coefficients.ci(model = lm.mod5, mod.vcov = lm.mod5$clus.vcov, baseline.var2 = baseline.male, baseline.int = baseline.female)
model6est <- sum.coefficients.ci(model = lm.mod6, mod.vcov = lm.mod6$clus.vcov, baseline.var2 = baseline.male, baseline.int = baseline.female)
modelGAMest <- sum.coefficients.ci(model = gam.mod.out, baseline.var2 = baseline.male, baseline.int = baseline.female)

female.ys <- rev(seq(0.5,3.5,0.5))
male.ys <- rev(seq(4,7,0.5))

male.col <- "gray"
female.col <- "black"
background.col <- "#F3F4E8"

for.plotting <- rbind(model1est, model2est, model3est, model4est, model5est, model6est, modelGAMest)
xlims <- range(c(for.plotting[,c("var1.lower","var1.upper","sum.upper","sum.lower")],-0.001))
lwd.choice <- 4

for.plotting <- for.plotting * 100
xlims <- range(c(for.plotting[,c("effect.size.var2.lower","effect.size.var2.upper","effect.size.int.upper","effect.size.int.lower")]))
lwd.choice <- 4
png("plots/responsiveness_effect_size.png",1200,550)
par(mar=c(4,4,0,11), xpd=F, cex=2, bg="transparent")
plot(for.plotting$effect.size.int, female.ys, xlim=xlims, ylim = c(0, 7.5), col= female.col, pch=19, bty="n", xlab="Effect size (%) relative to male minister baseline", yaxt="n", ylab="")
segments(x0 = for.plotting$effect.size.int.lower, x1 = for.plotting$effect.size.int.upper, y0= female.ys, col= female.col, lty=c(1:6,3), lwd=lwd.choice)
points(for.plotting$effect.size.var2, male.ys, col= male.col, pch=19)
segments(x0 = for.plotting$effect.size.var2.lower, x1 = for.plotting$effect.size.var2.upper, y0= male.ys, col= male.col, lty=c(1:6,3), lwd=lwd.choice)
abline(v=0, lty=3, lwd=lwd.choice)
par(xpd=T)
legend(x=max(xlims),y=6.5, legend=c("Male", "Female"), col=c(male.col, female.col), lty=1, bty="n", cex=1, lwd=lwd.choice)
legend(x=max(xlims),y=4, legend = c("Base model",expression(+ lambda[m0]),expression(+delta[t]), expression(+lambda["m0"]+delta[t]), expression(+lambda[m0] + delta[t] +lambda[m1]*t),expression(+lambda[m0] + delta[t] +lambda[m1]*t + lambda[m1]*t^2),"GAM"), lty=c(1:6,3), bty="n", cex=1, lwd=lwd.choice)
dev.off()

## Save the output of model 6 as a tex file to use in the write up

sink("latex/tables/responsiveness_effect_size.tex")
cat(paste0(round(for.plotting[6,]$effect.size.int),"\\% [",round(for.plotting[6,]$effect.size.int.lower),"\\%, ", round(for.plotting[6,]$effect.size.int.upper),"\\%]"))
sink()

### ################################################
### Debate-level models
### ################################################

formulas <- list("~ minister_gender*Gender",
                 "~ minister_gender*Gender + debate_department",
                 "~ minister_gender*Gender + as.factor(yearmon)",
                 "~ minister_gender*Gender + debate_department + as.factor(yearmon)",
                 "~ minister_gender*Gender + debate_department*yearmon_numeric + as.factor(yearmon)",
                 "~ minister_gender*Gender + debate_department*poly(yearmon_numeric,2) + as.factor(yearmon)")

out_debate <- model.data[,list(sim.words = mean(sim.words),
                 n = .N,
                 minister_gender = unique(minister_gender),
                 debate_department = unique(debate_department),
                 yearmon = unique(yearmon),
                 yearmon.numeric = unique(yearmon.numeric)),
           by = list(subsection_id, Gender)]

mod1 <- as.formula("sim.words ~ Gender*minister_gender")
mod2 <- as.formula("sim.words ~ Gender*minister_gender + debate_department")
mod3 <- as.formula("sim.words ~ Gender*minister_gender + as.factor(yearmon)")
mod4 <- as.formula("sim.words ~ Gender*minister_gender + debate_department + as.factor(yearmon)")
mod5 <- as.formula("sim.words ~ Gender*minister_gender + debate_department*yearmon.numeric + as.factor(yearmon)")
mod6 <- as.formula("sim.words ~ Gender*minister_gender + debate_department*poly(yearmon.numeric,2) + as.factor(yearmon)")

nreps <- 1000

lm.mod1.deb <- bootstrap.function(mod1, cluster = "debate_department", data = out_debate, weights = out_debate$n, run.boot = T, reps = nreps)
lm.mod2.deb <- bootstrap.function(mod2, cluster = "debate_department", data = out_debate, weights = out_debate$n, run.boot = T, reps = nreps)
lm.mod3.deb <- bootstrap.function(mod3, cluster = "debate_department", data = out_debate, weights = out_debate$n, run.boot = T, reps = nreps)
lm.mod4.deb <- bootstrap.function(mod4, cluster = "debate_department", data = out_debate, weights = out_debate$n, run.boot = T, reps = nreps)
lm.mod5.deb <- bootstrap.function(mod5, cluster = "debate_department", data = out_debate, weights = out_debate$n, run.boot = T, reps = nreps)
lm.mod6.deb <- bootstrap.function(mod6, cluster = "debate_department", data = out_debate, weights = out_debate$n, run.boot = T, reps = nreps)

out_debate$debate_department_factor <- as.factor(out_debate$debate_department)

gam.mod.out.deb <- gam(sim.words ~ Gender*minister_gender + s(yearmon.numeric, by=debate_department_factor) + debate_department_factor + as.factor(yearmon), data=out_debate, weights = n)
names(gam.mod.out.deb$coefficients)[1] <- "Constant"
names(gam.mod.out.deb$coefficients)[grep("GenderF:minister_genderF",names(gam.mod.out.deb$coefficients))] <- "Interaction"
gam.se.deb <- summary(gam.mod.out.deb)$se
names(gam.se.deb)[1] <- "Constant"
names(gam.se.deb)[grep("GenderF:minister_genderF",names(gam.se.deb))] <- "Interaction"

rename_coefs <- function(true.model){
  
names(true.model$coefficients)[1] <- "Constant"
names(true.model$se)[1] <- "Constant"
names(true.model$clustered.se)[1] <- "Constant"
names(true.model$boot.se)[1] <- "Constant"

names(true.model$coefficients)[grep("GenderF:minister_genderF",names(true.model$coefficients))] <- "Interaction"
names(true.model$se)[grep("GenderF:minister_genderF",names(true.model$se))] <- "Interaction"
names(true.model$clustered.se)[grep("GenderF:minister_genderF",names(true.model$clustered.se))] <- "Interaction"
names(true.model$boot.se)[grep("GenderF:minister_genderF",names(true.model$boot.se))] <- "Interaction"
return(true.model)
}

lm.mod1.deb <- rename_coefs(lm.mod1.deb)
lm.mod2.deb <- rename_coefs(lm.mod2.deb)
lm.mod3.deb <- rename_coefs(lm.mod3.deb)
lm.mod4.deb <- rename_coefs(lm.mod4.deb)
lm.mod5.deb <- rename_coefs(lm.mod5.deb)
lm.mod6.deb <- rename_coefs(lm.mod6.deb)

ministry.fe <- c("Ministry FEs","$\\times$","\\checkmark","$\\times$","\\checkmark","\\checkmark","\\checkmark","\\checkmark")
month.fe <- c("Month FEs","$\\times$","$\\times$","\\checkmark","\\checkmark","\\checkmark","\\checkmark","\\checkmark")
linear.tt <- c("Linear time trends","$\\times$","$\\times$","$\\times$","$\\times$","\\checkmark","\\checkmark","$\\times$")
quadratic.tt <- c("Quadratic time trends","$\\times$","$\\times$","$\\times$","$\\times$","$\\times$","\\checkmark","$\\times$")
gam.tt <- c("Flexible time trends","$\\times$","$\\times$","$\\times$","$\\times$","$\\times$","$\\times$","\\checkmark")

sink("latex/tables/ols_responsiveness_tables_debate_level.tex")
mod_stargazer(lm.mod1.deb, lm.mod2.deb, lm.mod3.deb, lm.mod4.deb, lm.mod5.deb, lm.mod6.deb, gam.mod.out.deb,
              keep=c("Gender","minister_gender","Interaction","Constant"), 
              order=c("Gender","minister_gender","Interaction","Constant"), 
              se=list(lm.mod1.deb$boot.se, lm.mod2.deb$boot.se, lm.mod3.deb$boot.se, lm.mod4.deb$boot.se, lm.mod5.deb$boot.se, lm.mod6.deb$boot.se, gam.se.deb), 
              covariate.labels=c("Female","Female Minister","Interaction","Constant"), 
              keep.stat=c("n","rsq","adj.rsq"), 
              dep.var.caption="", 
              dep.var.labels="\\emph{Responsiveness}", 
              add.lines = list(ministry.fe, month.fe, linear.tt, quadratic.tt,gam.tt),
              column.sep.width="0.25pt",
              label="tab:ols_cosine", 
              title="Responsiveness (OLS)", 
              no.space=T, model.names=F)
sink()


### ################################################
### Validation checks
### ################################################

# Are ministers more responsive to backbenchers than backbenchers are to ministers in QT debates? 

similarity <- speeches[grep("Oral Answers to Questions", speeches$parent)]

similarity$minister.responding <- !similarity$minister_in_debate & similarity$next.speech.minister

similarity$backbencher.responding <- similarity$minister_in_debate & !similarity$next.speech.minister

similarity$min.min <- similarity$minister_in_debate & similarity$next.speech.minister
similarity$back.back <- !similarity$minister_in_debate & !similarity$next.speech.minister 
min.mod3 <- lm(sim.words ~  minister.responding,data=similarity[back.back==F & min.min==F])

sink("latex/tables/minister_responsiveness_all_data.tex")
mod_stargazer(min.mod3, covariate.labels=c("Minister responding to backbencher"),label="tab:minister_response",dep.var.labels=c("\\emph{res}"),dep.var.caption="",title="Minister and backbench responsiveness",keep.stat=c("n","rsq"), no.space=T)
sink()
